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Abstract 

Metapopulation models provide the theoretical framework for describing disease spread between different populations 
connected by a network. In particular, these models are at the basis of most simulations of pandemic spread. They 
are usually studied at the mean-field level by neglecting fluctuations. Here we include fluctuations in the models 
by adopting fully stochastic descriptions of the corresponding processes. This level of description allows to address 
analytically, in the SIS and SIR cases, problems such as the existence and the calculation of an effective threshold for 
the spread of a disease at a global level. We show that the possibility of the spread at the global level is described 
in terms of (bond) percolation on the network. This mapping enables us to give an estimate (lower bound) for the 
pandemic threshold in the SIR case for all values of the model parameters and for all possible networks. 
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1. Introduction 

Modeling the spread of a disease between different 
populations connected by a transportation network is 
nowadays of crucial importance. Global spread of pan- 
demic influenza is an illustrative example of the impor- 
tance of this problem for our modern societies. In the 
context of pandemic spread and airline transportation 
networks, Rvachev and Longini [1] proposed a model, 
now called a metapopulation model, which was used to 
describe the spread of pandemic influenza [2, 3, 4, 5], 
SARS [6, 7], HIV [8], and very recently swine flu [9, 
10]. This model describes homogeneously mixed popu- 
lations connected by a transportation network. Despite 
its successes, very few theoretical studies are available. 
Numerical studies [11] showed that weight heterogene- 
ity played a major role in the predictability by creating 
epidemic pathways which compensate the various sce- 
nario possibilities due to hubs [7]. 

A major advance in our understanding of metapopu- 
lation models was done in the series of papers by Col- 
izza, Pastor-Satorras and Vespignani [12, 13, 14]. In 
particular, these authors studied the condition at which a 
disease which can spread in an isolated population (and 
has therefore a basic reproductive rate fio > 1 (see be- 
low)) can invade the whole network. In this respect, 
these authors defined an effective reproductive rate "R» 



at the global level and showed that its expression for a 
random network, in the SIR case with 7?o — > 1, and with 
a travel rate p, is under a mean-field approximation [12] 



K, = (Ko - 1) 



(k 2 ) - (k) pNa 



(1) 



where the brackets (•) denote the average over the con- 
figurations of the network, while a represents the frac- 
tion of infected individuals in a given city during the 
epidemic, 1/jU is the mean time an individual stays in- 
fected, and N is the mean population of a city. For 
< R„> \ the disease spreads over an extensive number of 
cities. For scale-free networks, the second moment (k 2 ) 
being very large, travel restrictions (i.e., decreasing p) 
have a very mild effect. Equation (1) relies on a number 
of assumptions: the network is uncorrelated and tree- 
like, 7?o is close to 1, and more importantly the number 
of infected individuals going from one city to another is 
assumed to be given by the constant d = pNa/p. This 
last expression however neglects different temporal ef- 
fects. Indeed, in the SIR case we can distinguish two 
different time regimes: the first regime of exponential 
growth corresponding to the outbreak, and the recov- 
ery regime where, after having reached a maximum, the 
number of infected individuals falls off exponentially. 
The number of infected individuals going from one city 
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to another thus depends on different factors and can be 
very different from the expression for d given above. 

In the present work, we consider the stochastic ver- 
sion of the metapopulation model. We investigate in de- 
tail the number of infected individuals going from one 
site to the other. This allows us to address the problem 
of the condition for propagation at a global scale with- 
out resorting to mean-field approximations or Kq — > 1 . 
We will show that the latter can be recast in terms of a 
bond percolation problem. 

So far most of the simulations on the metapopulation 
model were done in essentially two ways. The first ap- 
proach consists in using Langevin-type equations, with 
a noise term whose amplitude is proportional to the 
square root of the reaction term [15, 5]. In this approach, 
the evolution of the various populations (infected, sus- 
ceptibles, removed) is described by a set of finite differ- 
ence equations, with noise and with traveling random 
terms describing the random jumps between cities [5]. 
This level of description suffers from several technical 
difficulties when the finite-difference equations are it- 
erated: truncation of noise, non-integer number of in- 
dividuals, etc., which render its numerical implementa- 
tion difficult. Another possibility consists in adopting 
an agent-based approach, i.e., simulating the motion of 
each individual [5]. This simulation is free of ambigui- 
ties but requires a lot of memory and CPU time. 

In contrast, in the present work we will put forward an 
intermediate level of description, which is free from the 
technical difficulties of the Langevin-type approach (it 
does not require to iterate finite-difference equations), 
and which does not require large amounts of CPU. This 
approach consists in describing the process at the pop- 
ulation level, considering that individuals are indistin- 
guishable. In Statistical Mechanics, this level of de- 
scription is used for urn models or migration processes 
(see [16, 17] for reviews). Let us take for definiteness 
the example of the stochastic SIS epidemic process. For 
a single isolated city with population N, this process 
(see for example [18]) is described in terms of two vari- 
ables: S(t) (number of susceptible individuals) and /(f) 
(number of infected individuals), which evolve in con- 
tinuous time as 

(/ S) -» | (/ + h S ~ 1} Wi * ASI/N ' (2) 
[(7-1,5+1) with rate fil, 

where S (/) + I(t) = N is conserved, and with initial con- 
dition I(t = 0) = /o- The first reaction corresponds to 
the meeting of a susceptible individual with an infected 
individual, resulting in two infected ones. The second 
reaction corresponds to the recovery of an infected in- 



dividual into a susceptible one. The SIR case has an 
analogous definition (see Section 3). 

For two or more cities, we have to take into account 
the process of traveling between different cities. This 
diffusion process is described at the microscopic level 
by the jumping of individuals from one city to another. 
Let ;' and j be two neighboring cities. As described 
above, at the level of populations (i.e., such that the in- 
dividuals are indistinguishable), the process inside each 
city is given by (2), while the elementary traveling pro- 
cesses between these two cities are: 

(/,, /,) -» (7j - 1, 1 j + 1) with rate /?,/,■, 
(S h Sj) -» (St - l,Sj + 1) with rate p u Si, 1 ; 

where pij is the probability per unit time and per indi- 
vidual to jump from city i to city j. 

As said above, the numerical implementation of these 
processes is easy, and it allows us to make simula- 
tions with large numbers of large cities, in contrast 
with agent-based simulations. In the following we will 
restrict ourselves to uniform symmetric diffusion, de- 
fined by the constant rates pij — p if cities ; and j 
are neighbors. These rules are compatible with a sta- 
tionary state where all cities have the same population 
(/, + S i — Nj = N) on average. We will consider the 
initial condition /,■(* = 0) = 7o<5io, where Iq is the initial 
number of infected individuals, which all live in city at 
time t — 0. All the numerical simulations of the SIS and 
SIR cases will be performed with A = 0.3 and p = 0.1, 
so that ^o = 3 (see (6)). 

The outline of the paper is as follows. In Section 2 we 
consider the SIS case in different geometries, starting 
with the case of a single city, then proceeding to the case 
of two cities, and finally discussing the propagation on 
an infinite one-dimensional array of cities. As the num- 
ber of infected individuals in a given city saturates to a 
non-zero fraction of the total population, the metapopu- 
lation model will always experience a pandemic spread 
and ballistic propagation. It is however instructive to 
study this model for the following reasons: (i) it dis- 
plays the same phase of exponential growth as the SIR 
model, but with the great simplification that only two 
compartments (S and I) are present, allowing analyti- 
cal approaches; (ii) the SIS model has connections with 
the FKPP equation [19, 20], which is still currently the 
subject of intense activity. In Section 3 we consider 
the SIR case along the same lines. We successively 
consider the simple cases of one and two cities, before 
addressing the question of the propagation on various 
extended structures (one-dimensional chain, square lat- 
tice, Cayley tree, uncorrected scale-free network). The 
key difference with the SIS case is that the time integral 



2 



of the number of infected individuals in a given city is 
now finite. We demonstrate that the latter property im- 
plies the existence of an effective pandemic threshold 
for the spread of the disease over the whole network. 
We also obtain an estimate (lower bound) for this pan- 
demic threshold as a function of all the model parame- 
ters and of the bond percolation threshold on the under- 
lying geometrical structure. For scale-free networks our 
prediction is virtually identical to that of [12]. 

2. Metapopulation model in the SIS case 

In this section we discuss the metapopulation model 
in its stochastic form, on the example of the SIS pro- 
cess. We start with the case of a single city, then pro- 
ceed to the case of two cities, which will be the building 
block of our analysis of the propagation of the epidemic 
on extended structures. We finally discuss the ballistic 
propagation on a one-dimensional array of cities. 

2.1. Stochastic SIS model for one city 

We first present a discussion of the stochastic version 
of the SIS model for one city. 

Let N be the population of the city and /(f) the num- 
ber of infected individuals at time t. As can be seen 
from the defining reactions (2), the SIS process can be 
described by the single random variable 7(f), which ex- 
periences a biased continuous-time random walk on the 
interval (0,N), with variable rates which only depend 
on the instantaneous position 7(f) = k of the walker. 
This process is known as a birth-and-death process in 
the probabilistic literature [21]. Let us denote by Ak and 
Hk the jump rates of the walker, respectively to the right 
(k -» k + 1) and to the left (k -> k - 1). Here we have 



h = & 1 



Hk = ixk. 



(4) 



Hence the origin (k = 0) is absorbing while the right 
end (k = N) is reflecting. The walker moves in an effec- 
tive confining potential, with a restoring force bringing 
it back to the quasi-equilibrium position corresponding 
to the equality of the rates Ak and This position can 
be interpreted as the number of infected individuals at 
saturation. It reads 



/sal = N [ 1 - — I . (5) 

The so-called basic reproductive number is defined as 

(6) 



and we will hereafter focus on the case of interest 
where Hq is larger than 1, which is the condition for 
the occurrence of an outbreak at the level of a single 
city. The quasi-equilibrium state around 7 sat is however 
metastable, since the walker is deemed to be ultimately 
absorbed at the origin, which corresponds to the extinc- 
tion of the epidemic in the city. The characteristic ab- 
sorption time is however exponentially increasing with 
the population N of the city, as will be shown below. 

Let us denote by fk(t) = Prob(/(f) = k) the proba- 
bility for the number of infected individuals to be equal 
to the integer k at time t. This probability contains all 
the information on the one-time properties of the pro- 
cess (see [16, 17] for reviews). Its temporal evolution is 
given by the following master equations, characteristic 
of birth-and-death processes [22, 23]: 

dt 
d/o 
df 
dfy 
df 



: Ak-\fk-i + Hk+\fk+\ 



(A k +Hk)fk, 



(7) 



The nonlinear dependence of the rates on the position k 
renders their analytical study difficult. For example the 
evolution equation of the first moment reads 



(8) 



whereas the equation for {I 2 ) itself involves (7 3 ) and so 
on. However, eqs. (7) are easily implemented numeri- 
cally, for a given initial condition. In Figure 1 (top) we 
show a plot of (7) obtained by a numerical integration 
of eqs. (7). We also plot the same quantity as obtained 
from a numerical simulation of the SIS process. The 
data were obtained for A = 0.3, yu = 0.1, N = 100, and 
Io = 4, discarding histories leading to absorption. The 
two sets of data are indistinguishable one from the other. 
Finally we also show for comparison the deterministic 
(or mean-field) expression of 7(f) (see (10) below). Fig- 
ure 1 (bottom) shows five different realizations of the 
stochastic process for a single city, together with the de- 
terministic solution (10). 

The deterministic (or mean-field) approximation con- 
sists in neglecting any correlations. For example, in (8), 
it amounts to replacing (I 2 ) by (I) 2 . This approach is a 
priori legitimate if the population N of the city is large, 
so that the number of infected individuals 7(f) is also 
typically large enough to be considered as a determin- 
istic observable, instead of a random variable. This ap- 
proximation thus reduces the SIS model for one city to 
a deterministic dynamical system. Then (8) becomes 
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the following dynamical equation for the temporal evo- 
lution of 7(f): 



dl 



U-ti)l- 



dt r ' N 

the solution of which can be cast into the form 



1 

W) 



(9) 



(10) 



As demonstrated by Figure 1, the model exhibits 
an exponential growth regime followed by a saturation 
regime. The deterministic approximation gives an ac- 
curate global description of /(?), even for rather small 
populations. It however misses the effect of fluctuations 
in the process, that we now study in both regimes suc- 
cessively. 




• Numerical simulation 

— Numerical integration of the master equation 
■ ■■ Deterministic solution 
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Figure 1: (Color online). Top: Mean number of infected individuals 
in the stochastic SIS process for one city, conditioned on no extinction 
of the epidemic. Smooth red line: numerical integration of the master 
equation (7). Blue circles: Results of a numerical simulation. Dotted 
black line: deterministic solution (10). Bottom: Five different realiza- 
tions of the stochastic SIS process. Smooth black line: deterministic 
solution (10). (« = 3, N = 100, / = 4.) 



Exponential growth regime 

The regime of exponential growth formally corre- 
sponds to taking the infinite N limit at fixed time t. In 
this limit, (8) and (9) simplify, and their common solu- 
tion reads 

(I(t))=I ^ )t . (11) 



Furthermore, the rates become linear in k, i.e., Au - 
and iik - file, reducing the model to a solvable birth-and- 
death process [21, 23]. The master equations for the ft 
read 

^ = A(k - l)/*_! + fx(k + l)f k+l -(A + ii)kf k , 
dt 



d/o , 

The generating function 

G(s,f) = (s I ^) = J]s k f k (f) 



(12) 



(13) 



AMI 



satisfies 



(s - l)(As - fi)- 



dG 



(14) 



dG 

dt v " ~' v ~ *~' ds 

the solution of which is obtained with the method of 
characteristics and is 



^^l ^-D-^-^e-^ J ' (15) 

In particular, the probability for the number of in- 
fected individuals to be zero at time t, fo(f) = G(0, f), 
reads [23] 



Mt) = 



A -yue 



-{A-n)t 



(16) 



For t — * co we recover the well-known extinction prob- 
ability [24] 

1 



Pextinct = /()(°°) 



(17) 



If flo < 1 this extinction probability is equal to 1 . Hence 
the condition Hq > 1 ensures the possibility of a local 
outbreak at the level of a single city. The expressions of 
the moments of I(t) can be extracted from (15). We thus 
recover (11), while the variance reads 

ysrm= ^ii<m « /(0 > _7 0) . ( i8) 

In the late stages of the exponential growth regime, the 
relative variance is proportional to 1/Iq: 



var/ _ gp + 1 J_ 

W ~ %> - 1 v 



(19) 



Saturated regime 

After the phase of exponential growth is over, the 
mean number of infected individuals saturates at the 
value (1(f)) = 7 sat given by (5), as seen in Figure 1. 
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The time to reach saturation therefore scales as f sat « 

The saturated state is actually a metastable state, 
rather than a genuine stationary state. For finite N, the 
system will indeed always end in its absorbing state at 
k = 0. The distribution of the fluctuating number of 
infected individuals in the metastable state can be eval- 
uated by cutting the link from k — 1 to k — which is 
responsible for absorption. The process thus modified 
reaches an equilibrium state at long times. In the latter 
state the distribution fa satisfies detailed balance, that is 

Wfc+i/i+i = hfk, (20) 
where and p^ are defined in (4). This equation yields 



A 



1 (AO! (*o 



JVgCft/AO 



(21) 



fa tR (N-ky.\N 
where the large deviation function g(x = k/N) reads 

g(x) = x(ln^ - 1) - (1 - x)ln(l - x). (22) 

As expected, the distribution fa is peaked around k = 
Nx c = 7 sat for N large. Indeed the function g(x) takes its 
maximal value, 



= lnfto + - 1, 



(23) 



forx = x c - 1— l/9?o. We therefore have the exponential 
estimate fajfa ~ e Ngmia . The lifetime of the metastable 
state can then be estimated, in the spirit of the Arrhe- 
nius law, as r ~ 1 /fa . It is therefore predicted to grow 
exponentially with the population as 



tn 



Jlgm 



(24) 



For 7?o = 3 we have g max = 0.431945. 

The expression (21) also yields an estimate for the 
Gaussian fluctuations of 7(f) around its value 7 sat in the 
saturated state. Indeed, expanding g(x) to second order 
around x c , we obtain the estimate 



N 



(25) 



for the variance of the distribution of I in the quasi- 
stationary state, so that the reduced variance 



var/ <R<> 1 



</> 2 (Kb - I) 2 N 



(26) 



is proportional to 1 /N, whereas it was proportional to 
I Ik) in the growth phase (see (19)). These two results 



provide a quantitative confirmation that relative fluctua- 
tions around the deterministic theory become negligible 
in all regimes, as soon as the number of infected indi- 
viduals is large. Figure 2 shows a plot of var/, obtained 
by integration of the master equations (7). The data for 
small and large times are found to be in very good agree- 
ment with the estimates (18) and (25), respectively. 




Figure 2: (Color online). Variance of the number of infected indi- 
viduals in the stochastic SIS process for one city (Hq = 3, N = 100, 
Iq = 4). Continuous line: numerical integration of the master equa- 
tions (7). Red dashed line: prediction (18). Blue dot-dashed line: 
stationary value (25). 



2.2. Including travel: two cities 

Consider now two cities, numbered and 1 . The trav- 
eling process is described by (3). We assume symmetric 
diffusion: poi = pm = p, hence the system reaches a 
stationary state such that No = N\ — N. We choose the 
initial condition 7,(f = 0) = /o£;o- We are primarily in- 
terested in the time of occurrence of the outbreak of the 
epidemic in city 1 . This event is governed by the arrival 
of infected individuals traveling from city to city 1. 
Returns of infected individuals from city 1 to city can 
be neglected throughout. 

We start by analyzing the distribution of the first ar- 
rival time 1 1 of an infected individual in city 1 . Let us 
denote by Q{a, b) the probability that no infected indi- 
vidual exits from city in the time interval (a, b), with 
in particular Q(0, t) = Prob(?i > t). We have 



2(0, f + d?) = Q(0,t)Q(t,t + dt) 

= Q(0,t)(l-pI Wt), 



(27) 



where the expression in the last parentheses is the prob- 
ability that no infected individual exits from city in the 
infinitesimal time interval (f, t + dt). We thus obtain 



d6(0,f) 
dt 



- P I (t)Q(0,t), 



and so 



Prob(fi > t) = e 



-A(f) 



(28) 
(29) 
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with 



A(f) 



Jo 



dT/ (T). 



(30) 



An alternate derivation of this result Eqs. (29,30) is 
given in [25]. We thus conclude that the number N{t) 
of arrivals of infected individuals in city 1 in the time 
interval (0, f) is a Poisson process, for which the rate is 
itself a stochastic process, equal to dA/df = plo(t), as 
the integrated rate is A(f). Therefore 



Ait)" 

Prob(JV(r) = n) = e~ A(,) — — (n = 0, 1, . . .). 



(31) 



(h) 4 



1/pN 



■In 



X — fj, 



,A - (i pl 




A — fi 
N 



Pi 



A — ii 



* P 



Figure 3: Schematic (log-log) representation of the three different 
regimes for the mean arrival time (t[ ) as a function of p. 



This generalization of the Poisson process is known in 
the literature as the Cox process [26]. Hereafter, in order 
to simplify the presentation, we will implicitly assume 
that probabilities are conditioned on a single stochastic 
history Io(t). 

The mean first arrival time can easily be derived 
from (29). It has the simple form 



(h) 



r 

Jo 



dfe- A(,) . 



(32) 



This exact expression can not be written in closed form, 
even if 7o(f) is given its deterministic value (10). The 
scaling properties of (t\ ) at large N can however be ob- 
tained using the following simple arguments. In the 
large p regime, we expect {t\) ~ l/(p/o), as the number 
of infected individuals has hardly changed from its ini- 
tial value Iq in time {t\ ). In the opposite regime, where p 
is very small, we have Io(t) ~ ^sat (see (5)), and thus 
(t[) as l/(p/ sat ) ~ l/(pN). In the intermediate regime, 
city is in its phase of exponential growth, and the mean 
arrival time can be estimated by imposing that A((fi )) is 
of order unity. We thus have 



<fi> 



1 



In 



A - pi 



A- p pl 



(33) 



The crossovers between these three regimes occur at 
p\ « (A - p)/Io and p2 ~ (A - p)/N. Figure 3 sum- 
marizes the above discussion, whereas actual data will 
be shown in Figure 5. 

We now turn to the distribution of the outbreak 
time fj, defined as the time of the arrival in city 1 of 
the first infected individual who induces an outbreak in 
city 1 . If the travel rate p is small enough, events where 
two or more infected individuals coming from city 
are simultaneously present in city 1 can be neglected. 
Each infected individual entering city 1 has therefore a 
chance l/fto (see (17)) to disappear before it induces 
an outbreak. On the other hand, the distribution of the 



arrival time ti of the k-th infected individual reads 



/t-i 



Probfo > = Prob(JV(f) <k) = e~ A(,) V — — . (34) 



n=0 



The probability that no outbreak occurred up to time t is 
thus given by 

Prob(f; > t) = Y(l - l/Ko)ftJ-*Probfo > 



k>\ 

-(l-l/R )A(f) 



(35) 



Comparing this expression with the corresponding one 
for the first arrival (29) reveals that taking into account 
extinction amounts to renormalizing p to an effective 
rate p*, and accordingly A(f) to A*(f), multiplying both 
quantities by the scaling factor 1 - l/'Ro- 

p * = [ l -k) p > = W (36) 

Within this approach, the distribution of the outbreak 
time in city 1 reads 



Prob(f[ > f) = e 



-A*(0 



(37) 



2.3. Ballistic propagation on a one- dimensional array 

We now consider the stochastic SIS model in the one- 
dimensional geometry of an infinite array of cities with 
initial populations N. We still choose a symmetric travel 
rate p, hence, in the course of time, these cities remain 
equally populated on average. 

The model is observed to reach very quickly a ballis- 
tic propagation regime, where the epidemic invades the 
whole array by propagating a ballistic front moving at a 
well-defined finite velocity V. Such a front is illustrated 
by the space-time plot of Figure 4. 

The velocity V of the front can be estimated as fol- 
lows. The time T„ at which the epidemic reaches city n 
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Figure 4: Propagation of a ballistic front in the SIS case (Ko = 3, 
p = 0.01, N = 100) shown in the usual space-time representation 
where the y-axis represents the time (in arbitrary units) and where the 
jc-axis represents the index of cities. 



can be recast as 



T n = Y j f 1 (i), 
1=1 



(38) 



where f*(;) is the outbreak time in city i, i.e., the ar- 
rival time of the infected individual which will trigger 
the outbreak in that city, with the origin of times being 
set to the outbreak time in the previous city i— 1. Model- 
ing the propagation from city i— 1 to city ; by the case of 
two cities studied above, we thus predict that the inverse 
velocity l/V = lim(r„/n) reads 



1 

V 



<0, 



(39) 



where (?*) is the mean outbreak time of the problem of 
two cities, with the natural initial condition Io — 1 . 

In Figure 5 we show a comparison between data for 
the reciprocal of the ballistic velocity V and for (fj), the 
mean outbreak time for two cities, with Iq — 1. The 
latter was measured in a Monte Carlo simulation as (t i ) 
with p renormalized to p* (red symbols), and computed 
using the analytic prediction (32), where 7o(f) has the 
deterministic expression (10), again with p renormal- 
ized (blue curve). The agreement validates the approx- 
imation (39). The prediction (43), (44) of the discrete 
FKPP theory is also shown. 

It is indeed worth confronting the above analysis with 
yet another approach, based on the analysis of the fol- 
lowing deterministic equation for the densities of in- 
fected individuals p,(f) = Ii(t)/N: 



dpi(t) 
dt 



= Api(\ - pi) - ppi + p{pi+\ + pi-\ - 2pi). (40) 



10" 

p 



Figure 5: (Color online). Upper black curve: Reciprocal of the ballis- 
tic velocity V against p (Rq = 3, N = 100, p = 0.01). Blue curve with 
red symbols: (t[) of the two-city problem (analytic prediction (32)) 
and numerical simulation, with Iq = 1 and p renormalized. Lower 
smooth green curve: prediction (43), (44) of FKPP theory. 



it follows that the function fix) obeys the differential- 
difference equation 

-Vf(x) = Af(x)(l-f(x))-p.f(x) 

+ p(f(x+\) + f(x-\)-2f(x)). (41) 

If one assumes that the solution decays exponentially as 
fix) ~ e~°~ v for x — * co, the velocity is found to depend 
continuously on cr, according to 



V 



A — p + 2p(cosh cr — I) 



(42) 



For a localized initial condition, the actual velocity of 
the front is known [27, 28] to be obtained by minimizing 
the above expression with respect to the spatial decay 
rate cr. This velocity reads 



V 



(A - fi) sinh cr 



cr sinh cr + 1 - cosh cr 
and it is reached for 

A-p 



2(cr sinh cr + 1 - cosh cr) 



(43) 



(44) 



Both above equations give a parametric expression for 
the ballistic velocity V. The usual prediction of the con- 
tinuum FKPP equation, i.e., V = 2{{A - p)p) 1 ^ 2 , is re- 
covered as cr «c 1, i.e., p » A-p. The regime of 
most interest in the present context is the opposite one 
(cr » 1, i.e., p <K A - p), where we have 



1 

V 



J- tail. 
A-p pe 



(45) 



This equation is the discrete version of the FKPP equa- 
tion [19, 20]. Looking for a traveling- wave solution of 
the form p,(f) = f(i - Vt), propagating with velocity V, 



This regime is in correspondence with the estimate (33) 
which holds in the intermediate growth regime. Fig- 
ure 5 shows that the discrete FKPP theory provides a 



good, albeit not quantitative, description of the overall 
behavior of the ballistic velocity. 

To conclude, in the one-dimensional SIS model, the 
epidemic always spreads ballistically, with a non-zero 
velocity V. We expect this result to hold on any ex- 
tended structure, since it is merely a consequence of the 
existence of a very long-lived saturated regime where a 
finite fraction of individuals are infected. In a metapop- 
ulation model, this will for sure trigger an infection 
which will spread over the whole network as long as 



3. Metapopulation model in the SIR case 

We now turn to the case of the SIR model for the 
growth of an epidemic in a given city. In contrast with 
the SIS model, the number of infected individuals in a 
given city now falls off to zero at large times. As a con- 
sequence, if the traveling rate is too small compared to 
the typical inverse duration of the epidemic phase, the 
epidemic may die out in a city before propagating to the 
neighboring ones. This phenomenon will lead to the ex- 
istence of a non-trivial threshold /? t h for the travel rate, 
below which the disease will not spread in the system. 
We will show that the problem can be recast in terms of 
a bond percolation problem. This picture greatly simpli- 
fies the analysis, and leads to an estimate (lower bound) 
of the pandemic threshold. 



3.1. Stochastic SIR model for one city 

We first investigate the stochastic SIR process in the 
case of one isolated city of population N. As in the 
SIS case, the process is described at the population level 
(i.e., the individuals are indistinguishable), but we now 
have three variables S(t), /(?), and R(t) (number of re- 
covered individuals), with S (t) + 1(f) + R(t) = N. The 
stochastic process is described by the reactions 



The deterministic equations describing the SIR pro- 
cess read 



(S,I,R) 



\(S - 1,7 + l,R) with rate AS I/N, 
(S,7- 1,7?+ 1) with rate //7, 



(46) 



with initial condition 1(0) = 7o and R(Q) = 0. The sec- 
ond reaction corresponds to the recovery of an infected 
individual, that stays immune later on. The master equa- 
tions for the process can be written in terms of two in- 
dependent integer random variables, say S and I, and 
integrated numerically. 



dS_ 
~dt 
dl 



A- 



.57 

' N' 



SI 

A ul, 

dt N 

dR 

— = ul. 

dt H 



(47) 



At variance with the case of SIS, the above equations 
cannot be solved in closed form. Figure 6 shows a com- 
parison between the numerical solution of these equa- 
tions and a numerical simulation of the stochastic pro- 
cess. The upper panel shows five different trajectories 
in the I-R plane. We observe relatively important fluc- 
tuations around the deterministic solution. In particu- 
lar the epidemic stops after a finite random time. The 
lower panel demonstrates that mean quantities are how- 
ever well described by the deterministic approach. 




Figure 6: (Color online). SIR model for one city. Top: five differ- 
ent realizations of the process plotted in the I-R plane. Smooth black 
line: solution of the deterministic equations (47) (N = 1000). Bottom: 
densities of susceptible (black), infected (red), and removed (green) in- 
dividuals as a function of time. Numerically measured mean densities 
(full) and corresponding deterministic solutions (dashed) are hardly 
distinguishable (%> = 3, N = 100, 7 = 10). 

A quantity of central interest for the sequel is the fi- 
nal size of the epidemic, defined as the number R(oo) of 
recovered individuals after the epidemic has stopped: 



R(oo) 



Jo 



dr/(r). 



(48) 



We refer the reader to [29, 30, 3 1 ] for studies on the dis- 
tribution of this quantity. For a city with a large popula- 
tion (N » 1), R(oa) can be approximately determined in 
the framework of the deterministic approach. It follows 
from (47) that 

^ = -*oi (49) 



dR 

from which we obtain 



N 



S(oo) R(oo) 

In = -%) ■ 



S{0) 



N 



(50) 



Using the initial condition given above, and the final 
condition 7(°°) = 0, i.e., S (°°) + R(oo) = N, we obtain a 
transcendental equation relating the densities z'o = Iq/N 
and roo = R(oo)/N: 



l-r 00 = (l-/o)e^ ''", 
which reduces for 7(0) <K N, i.e., io «: 1, to 

The density r^, starts rising linearly as 

roo * 2(7? - 1) Wo -» 1), 



(51) 



(52) 



(53) 



and reaches unity exponentially fast for 7?o — > oo, as 
^ w 1 - e~^°. For 7?o = 3, the value used in numerical 
simulations, we have r M = 0.940479. 

3.2. Including travel: two cities 

As in the SIS case, it is natural to start the study of 
propagation by the case of two neighboring cities. The 
SIR process inside each city is described in terms of 
three random variables Si, and R, (i = 0, 1), with 
Si + Ij + Ri = Ni, which evolve stochastically as in (46). 
The traveling process between these two cities is given 
by the reactions 

(Xi, Xj) -> (Xi -l,Xj + l) with rate ptjXi, (54) 

where X stands for S , / or R. We again assume symmet- 
ric diffusion. A typical history of the system is shown 
in Figure 7. 

The first arrival time t\ of an infected individual in 
city 1 is still distributed according to (29) and (30). 
There is however a crucial difference between SIS and 
SIR. As already discussed, in the present case of SIR, 
the integrated rate A(f) converges to a finite limit A^ 
as t — > oo. For a large population N, the deterministic 
approach yields 

Npr m 



Aoo = 



(55) 




Figure 7: (Color online). A typical history of the SIR process with 
travel between two cities: plot of the densities of susceptible, infected, 
and removed individuals versus time for cities and 1 . Color code as 
in Figure 6. 



As a consequence, there is now a non-zero probabil- 
ity exp(-Aoo) that no infected individual travels from 
city to city 1 during the whole epidemic in city 0. In 
other words, the event that at least one infected individ- 
ual reaches city 1, i.e., that the time t\ is finite, only 
occurs with probability 



n = 1 



1 



(56) 



The distribution of t\ is said to be defective [21]. 

Taking into account the fact that single infected in- 
dividuals only trigger an outbreak with probability 1 - 
1 fRo, the distribution of the outbreak time f* in city 1 is 
still given by (37), as a result of the renormalization pro- 
cedure which led us to replace p by the effective rate p* 
(see (36)). In particular, the probability of occurrence 
of an outbreak in city 1 reads 



IT = 1 



1 



(57) 



This probability is non-trivial, i.e., less than one, in con- 
trast with the SIS case. 

3.3. Propagation on extended structures 

We can now consider the general case where the cities 
are connected so as to form a network or any other kind 
of extended structure (one-dimensional array, regular 
finite-dimensional lattice, regular tree). Individuals can 
travel by performing diffusion along the links of the net- 
work, allowing thus the disease to spread over different 
cities. 

The probability IT given in (57) can be interpreted as 
the probability that the disease propagates through one 
given link from a city to one of its neighbors. In the SIS 
case, since the integrated rate A(f) diverges with time, 
the probability II* is trivially equal to unity. In contrast, 
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for the SIR process, this probability is less than unity, 
so that the disease can stop invading the network. 

For a given network, denoting the bond percolation 
threshold by p c , the condition for the disease to spread 
is therefore 

n* > p c . (58) 

This equation defines the pandemic threshold: for N 
large, (58) yields p > p t h, where 



Pth = 



A*|ln(l-j> c )| 
Ml - l/Ko)r« 



(59) 



This static prediction is not claimed to give an exact 
value for the pandemic threshold. It can however be ar- 
gued to provide a lower bound for the threshold, which 
is also meant as a reasonable and useful estimate. In- 
deed the picture of static percolation, where links are 
occupied with constant probability II*, corresponds to 
an ideal infinitely slow propagation, where the epidemic 
can take an infinitely long time to cross some of the 
links. In real situations, propagation is rather observed 
to take place with finite velocity. This velocity can be 
thought of as providing a time cutoff T for propagation 
across every single link of the network. The effect of 
this cutoff time is to decrease the integrated rate from 
A M to the smaller value A(T), and hence to increase the 
threshold value of p from the theoretical prediction (59) 
to a higher value. 

Let us now illustrate this result on some examples of 
networks and other extended structures. 




Figure 8: (Color online). Fraction of infected cities at a very large 
time versus travel probability p on the Cayley tree (k = 3, ; = 5 X 10 5 , 
Rq = 3, N = 100). The red dashed line points toward a threshold 
value near p t h ss 2.1 X 10~ 3 . 



On a 50x50 array, the data of Figure 9 clearly demon- 
strates the existence of a threshold near pth ~ 1 ■ 1 x 10" 3 , 
in good quantitative agreement (within ten percent, say) 
with the predicted static value. In figure 10 we provide 
an example of the observed shape of the infected region 
for a large but finite time. 




0.002 0.004 
P 



Cayley tree 

We consider first the geometry of a regular Cayley 
tree (or Bethe lattice) with coordination number k, so 
that p c = l/(k- 1). 

For the simplest example of k — 3, and an epidemic 
initially located at the center of a tree with 15 genera- 
tions, Figure 8 shows a plot of the fraction of infected 
cities at a very long time. We indeed observe a pan- 
demic threshold near pth ~ 2.1 x 10~ 3 . For p c = 1/2, 
H = 0.1, <Ro = 3, hence r^ = 0.940479, and N = 100, 
the predicted static threshold is p t h = 1.105 x 10~ 3 . The 
observed threshold is thus some two times larger than 
the static one. 

Square lattice 

We now consider the two-dimensional situation of 
propagation on the square lattice. This lattice also has 
bond percolation threshold p c = 1 /2, so that the above 
static threshold of p th = 1.105 x 10~ 3 still holds for the 
same parameter values. 



Figure 9: (Color online). Fraction of infected cities at a very large time 
versus travel probability p on the square lattice. Same parameters as 
in Figure 8. The red dashed line points toward a threshold value near 
p th * 1.1 xlO- 3 . 



One -dimensional array 

We finally address the case of propagation along an 
infinite array of cities. The percolation threshold in the 
one-dimensional case is p c — 1. Our analysis therefore 
suggests that the disease will always stop, after having 
invaded only a finite range of typical size but not 
the whole array. The static percolation approach sug- 
gests that £ diverges as IT — > 1 in the same way as the 
static correlation length in percolation theory, namely 
^ ~ 1/(1 - IT). We thus find the exponential growth 



£~exp 



P N(l - l/floK, 
A* 



(60) 



We indeed observe a symmetric ballistic front on both 
sides of the seed. Each branch suddenly dies at a ran- 
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Figure 12: (Color online). Cumulative distribution of the location m 
of the rightmost infected city (N = 100). The slope of the fitted red 
line yields £; = 12.4. 



Figure 10: Typical infected region on the square lattice after a large 
but finite time (system size 50 X 50, N = 100, p = 0.01 » pth). 



dom time, and hence in a randomly located city, say 
number m. Such a front is shown in Figure 1 1 . For given 
parameter values, the right stopping point m is observed 
to be exponentially distributed, with a cumulative dis- 
tribution falling off as exp(-m/£) (see Figure 12). This 
exponential distribution confirms our intuition that a lo- 
cal mechanism is responsible for the arrest of the prop- 
agation. The p and N dependence of the characteristic 
length £ thus measured is shown in Figure 13. The data 
demonstrate an increase of £ with both p and N. The 
agreement with the static prediction however remains 
very qualitative. 



Figure 11: (Color online). Propagation of a ballistic front in the SIR 
case (Hq = 3, p = 0.01, N = 100). Time (arbitrary units) versus index 
of the city. The front is observed to stop at a random time. 



Scale-Free networks 

We close up our analysis of the propagation on a 
SIR epidemic by considering the case of a scale-free 
network. For an uncorrelated network, the percolation 
threshold is known to be given by p c = (k)/((k 2 ) - (k)) 
(see (A. 7)). As most scale-free networks are such that 
(k 2 ) » (k), we have p c ss (k}/(k 2 ) « 1, so that our 




Figure 13: (Color online). Characteristic length f versus population 
N of one city. Colors correspond to various values of p. 



static prediction (59) for the pandemic threshold reads 
H (k) 1 



Pth 



N(k 2 )(l -l/KoK 



(61) 



In the same regime, it is worth recasting the result 
of [12], where the expression for 7?» is recalled in (1), 
in the form of a threshold value p± of the travel rate. 
Using the fact that a - r^ for SIR, we obtain 



Pth 



{k) 



Pth- 



(62) 



Our result thus only differs from that of [12] by an 
inessential multiplicative factor involving the mean de- 
gree (k) and the reproductive rate The coincidence 
is all the more striking that both approaches are largely 
different. Both expressions agree to predict that the 
threshold diverges quadratically as 1/(^0 - l) 2 as the 
basic reproductive rate 7?o — > 1. Indeed r m vanishes 
linearly (see (53)). 

4. Conclusion 

In this paper we have considered the metapopulation 
models for the SIS and SIR processes in their stochas- 
tic versions. We have put forward a description con- 
sidering individuals as indistinguishable, but keeping 
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the full stochastic character of the populations in the 
various compartments. This intermediate level of de- 
scription, along the lines of the theory of urn models 
and migration processes, allows very efficient numer- 
ical simulations. We have considered with great care 
temporal effects in the spread of a disease from one city 
to another and investigated its propagation at a global 
scale on scale-free networks and other extended struc- 
tures such as regular trees and lattices. 

For the SIS case we always observe a ballistic prop- 
agation. For the SIR case, even after an infinitely long 
time, the disease only gets transmitted with probabil- 
ity II* < 1 through every link of the network. This 
picture allows us to rephrase the possibility of global 
invasion as a static bond percolation problem. The re- 
sulting prediction is an estimate (lower bound) for the 
pandemic threshold, expressed as a threshold value p t h 
for the travel rate p. In the case of scale-free networks, 
our prediction is virtually identical to that of [12]. Our 
approach is however not limited to the case Rq — » 1, and 
takes into account both temporal and topological fluctu- 
ations. For other geometries (Cayley tree, square lat- 
tice), our static prediction yields a reasonable estimate 
for the pandemic thresholds observed in numerical sim- 
ulations. 

As in most previous theoretical studies, we assumed 
that the populations of all cities were equal and that the 
traveling rate per link and per individual was constant. 
It would be interesting to extend our result and to test for 
the relevance of travel and/or population heterogeneities 
on the existence of a pandemic threshold. 

Finally, it is worth putting the present work in a 
broader perspective. The relationship between epidemic 
spreading and percolation has already been discussed in 
the context of contact networks. Table 1 summarizes the 
parallel between the contact network and the metapop- 
ulation approaches. It has been emphasized first by 
Grassberger that the spreading of an epidemic over a 
contact network can be mapped onto directed percola- 
tion [32]. It was then realized that, as a general rule, 
epidemic models without immunization belong to the 
universality class of directed percolation, whereas those 
with immunization, such as SIR, belong to the univer- 
sality class of dynamical percolation (see [33] for a re- 
view). In the latter case, clusters of immune individu- 
als have the same critical properties as usual percolation 
clusters. The effective description of epidemic spread- 
ing put forward in this work demonstrates that epidemic 
spreading is also intimately related to static percolation 
at the metapopulation level. 



Network 


Nodes 


Links 


contact 


individuals 


contacts 


metapopulation 


cities 


travel 



Table 1 : A comparison between models of epidemic spreading defined 
on contact networks and at the metapopulation level. 



Appendix A. Pandemic threshold on a network 
with arbitrary travel rates 

In this appendix we determine the pandemic thresh- 
old p t h on an uncorrected network where travel is de- 
scribed by arbitrary hopping rates. 

We consider an uncorrected network, denote by Pk 
be the degree distribution of the nodes, and use the 
framework of degree classes, assuming that all nodes 
with given degree k are equivalent. Travel is defined by 
the rate pu for an individual to hop from a node with de- 
gree A: to a neighboring node of degree I. The stationary 
populations Nk obey the balance equation 

N^PiPki = J]PiPikNi, (A.l) 



where 



kP k 



(A.2) 



is the degree distribution of a neighboring node of 
a given node. For instance, uniform rates pu = p 
yield uniform populations Nk = N, whereas the rates 
Pu = p/k of ordinary random walk yield populations 
Nk = kN. More generally, separable (i.e., factorized) 
travel rates of the form pu = akbi yield stationary pop- 
ulations Nk = N bk/ak- 

The basic quantity in the metapopulation model is the 
probability that the epidemic will propagate (in an in- 
finitely long time) from a node of degree A: to a node of 
degree I. Dropping the star for conciseness, this proba- 
bility reads 



Uki = 1 - exp - 



N kP ki(l - 1/KoK 



P 



(A3) 



In order to determine the static pandemic threshold of 
SIR on the network, we are thus led to consider the 
problem of directed bond percolation on the network 
defined by the above probabilities Ylki that an oriented 
link is open from a node of degree A: to a node of de- 
gree I. Let us introduce the probabilities Qu that an ori- 
ented link from a node of degree k to a node of degree I 
leads to a node belonging to the giant component. These 
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quantities obey 

Qu = tiki 



where v runs over the (/ - 1) neighbors of the node of 
degree I which are not the initial one and m(v) are their 
degrees. Near the threshold, the probabilities Q k i are 
expected to be small. Linearizing the above relation, we 
obtain 

Qki = (l-l)TlkiYjP m Qi m . (A.5) 

m 

Setting Qu = Yl^iGi/Pi, the above condition reads G/ = 
Jj m M bn G m , with 



k(k - 1) 

(k) 



(A.6) 



The percolation threshold is therefore given by the con- 
dition that the largest (Perron-Frobenius) eigenvalue of 
the positive matrix M equals unity. 

The usual percolation problem corresponds to the 
case where every link is occupied with probability p, 
i.e., Tla = P for all k and I. We thus recover the known 
result 

* - why (A ' 7) 

In particular, for a regular Cayley tree with coordination 
number k, we have 



1 



1 



(A.8) 



The reasoning leading to (A.4) can be extended to other 
cases. Let us mention the example of a bipartite tree, 
where nodes with degree k\ are neighbors of nodes of 
degree ki- We then obtain 



1 



vc*i - ix*2 - 1) 



(A.9) 



For separable probabilities, of the form ILj = q k r\, we 
obtain the following threshold condition 



(k(k - \)q k r k ) 
<*> 



= 1. 



(A. 10) 



In a more general case, we have to resort to a numerical 
computation of the largest eigenvalue of the matrix M. 
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